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Abstract 

Evolutionary theory is sufficiently well developed to allow for short-term prediction of 
evolutionary trajectories. In addition to the presence of heritable variation, prediction 
requires knowledge of the form of natural selection on relevant traits. While many studies 
estimate the form of natural selection, few examine the degree to which traits evolve in 
the predicted direction. In this study we examine the form of natural selection imposed 
by mantid predation on wing size and shape in the fruitfly, Drosophila melanogaster. 
We then evolve populations of D. melanogaster under predation pressure, and examine 
the extent to which wing size and shape have responded in the predicted direction. We 
demonstrate that wing form partially evolves along the predicted vector from selection, 
more so than for control lineages. Furthermore, we re-examined phenotypic selection after 
~30 generations of experimental evolution. We observed that the magnitude of selection 
on wing size and shape was diminished in populations evolving with mantid predators, 
while the direction of the selection vector differed from that of the ancestral population 
for shape. We discuss these findings in the context of the predictability of evolutionary 
responses, and the need for fully multivariate approaches. 

Keywords: Predation, natural selection, experimental evolution, Drosophila, wing shape, 
geometric morphometries, adaptation, multivariate selection. 
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Introduction 

22 Biologists measure natural selection to help identify agents of selection, to infer how 

current phenotypes were influenced by past selection and to predict future evolutionary 
24 outcomes. Since the publication of the seminal work by Lande and Arnold Lande and 

Arnold (1983), considerable effort has gone into measuring the form, magnitude and vari- 
26 ability of phenotypic selection (Kingsolver et al., 2001; Hoekstra et a!., 2001; Kingsolver 

and Diamond, 2011; Siepielski et al., 2009; Morrissey and Hadfield, 2012). However, ad- 
28 ditional factors influence the trajectory of evolution such as correlational selection, the 

stability of the selective function (Brodie III, 1992), as well as the genetic architecture of 
30 the traits (Hansen and Houle, 2008; Agrawal and Stinchcombe, 2009; Kirkpatrick, 2009), 

making such predictions difficult in practice. In this study we address this predictability, 
32 by investigating the extent to which experimental populations of Drosophila melanogaster 

subject to predation risk, evolve along the trajectory predicted from phenotypic selection. 
34 Studies have investigated how closely populations evolve along the direction predicted 

from the multivariate breeders equation, (as a function of both selection and the ge- 
36 netic variance-covariance matrix) (Schluter, 1996; Hansen and Houle, 2008; Agrawal and 

Stinchcombe, 2009; Higgie and Blows, 2008; Hunt et al, 2007; Mcguigan et al, 2005; 
38 Walsh and Blows, 2009; Blows et al., 2004; Simonsen and Stinchcombe, 2010). Yet in 

only a handful of cases has selection been observed for a sufficient amount of time (be- 
40 yond a few generations) to evaluate these predictions. Furthermore, the ecology and 

natural history of many organisms limits us to estimating phenotypic selection, generally 
42 over just a few generations (but see Ozgul et al., 2009; Grant and Grant, 2002, 2006). 

For some selective agents like predation, the organism is consumed, prohibiting (at least 
44 in the field) the measurement of many traits that are targets of selection. As a result, 
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it may be challenging to predict the evolutionary trajectory of some traits involved with 
anti-predator activity. This might suggest a pessimistic view of our ability to predict the 
selective response in natural systems. 

Despite these issues, convergent and parallel evolution are often observed among pop- 
ulations, suggesting that persistent and predictable selection may be relatively common 
(Conte et al., 2012), even if it is difficult to measure. While estimates of the strength of 
viability selection suggest it may be weaker than for other fitness components (Hoekstra 
et al., 2001; Lind and Cresswell, 2005; Ajie et al., 2007), repeated evolution of simi- 
lar morphologies in response to predation for several fish species (O'Steen et al., 2002; 
Langerhans et al., 2004; Dayton et al., 2005; Langerhans and Makowicz, 2009) suggests 
a strong and consistent regime of selection. Similar results have also been observed for 
shell morphology among populations of snails in apparent response to predation (Auld 
and Relyea, 2011; DeWitt et al., 2000, 1999). When selection is relaxed by the removal 
of predators, even for just a few generations, trait means have been shown to change 
dramatically (Reznick et al., 1990, 1997; Reznick and Ghalambor, 2005), consistent with 
predation maintaining trait values in the face of potentially antagonistic selective effects. 
The prevalence of diverse, and often costly, traits that mediate interactions with predators 
suggests that predation profoundly influences fitness. 

In this study we investigate how multivariate wing form of Drosophila melanogaster 
evolves along the trajectory predicted by initial estimates of phenotypic selection in re- 
sponse to predation by mantid nymphs (Tenodera aridifolia sinensis). This novel ex- 
perimental system has a rather rare (but see Svensson and Friberg, 2007; Kuchta and 
Svensson, 2014) and useful attribute in which the wings are not consumed when the fly 
is captured by its mantid predator (figure 1A). This allows us to collect the wings from 
both surviving and consumed flies to estimate the form and magnitude of natural selection 
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on both size and shape. Multivariate shape provides a robust framework for evaluating 
evolutionary predictions. While size only varies along one axis, a high dimensional rep- 
resentation of shape is less likely to change in the predicted direction by chance alone. 
This enables us to make clear quantitative comparisons of the degree of similarity between 
predicted and observed response to selection (Pitchers et al., 2013). It also extends a well 
developed genetic system for use in studies of predator-prey interactions. 

Wing size and shape in Drosophila have been used as a model system for evolution 
(Gilchrist and Huey, 2004; Gilchrist et al., 2004; Huey et al., 2000; Gilchrist and Par- 
tridge, 1999; Weber, 19906; Mezey and Houle, 2005; Pitchers et al., 2013), genetics, and 
development (Dworkin and Gibson, 2006; Houle and Fierst, 2013; Palsson and Gibson, 
2000). There is substantial segregating variation for wing size and shape, with some vari- 
ants mapped (Weber et al., 1999; Palsson, 2004; Zimmerman et al., 2000; Mezey et al., 
2005; Mckechnie et al., 2010; Dworkin et al, 2005; Palsson et al., 2005). Studies have 
demonstrated that genetic variation is available along many dimensions of wing shape 
(Mezey et al., 2005). Using artificial selection, it has been demonstrated that this varia- 
tion can be selected upon (Houle et al., 2003; Weber, 19906, 1992, 1990a; Palenzona and 
Alicchio, 1973; Rochetta and Palenzona, 1975). Yet little is known about the selective 
agents influencing variation for wing form (but see (Hoffmann et al., 2007; Menezes et al., 
2013)) or the potential functional role it plays in avoiding predation. 

We quantified the magnitude and direction of selection on wing size and shape in an 
outbred population. We then allowed replicates derived from this population to experi- 
mentally evolve under episodic selection with mantids or under predator-free conditions. 
For the evolved populations we quantified changes in wing form, with particular focus 
on the direction of change, relating it to the vectors of phenotypic selection predicted 
from the base population. We demonstrate that while evolution of wing shape for the 
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predator populations is more aligned with the initial vector of selection than are the con- 
96 trols, not all of the change is in the direction predicted by the initial vector of selection. 

Despite observing consistent directional selection on wing size, we observed considerable 
98 divergence in its evolutionary response. Furthermore, we measured phenotypic selection 

on the evolved populations, and demonstrate that the magnitude of selection on both 
ioo size and shape is substantially diminished in the populations exposed to predation and is 

distinct from the initial vector of predicted selection. We discuss these results within the 
102 context of how populations change along a fitness surface, the importance of unmeasured 

traits and the degree of repeatability to agents of selection. 

104 Materials and Methods 

Base Populations 

we We used an advanced intercross with 100 inbred lines to generate a synthetic outbred 

ancestral population referred to as the base. The inbred lines were derived from two pop- 
loe ulations of wild Drosophila melanogaster collected in fruit orchards in Maine and North 

Carolina (Goering et al., 2009; Reed et a!., 2010). Flies were round robin intercrossed 
no for three generations and then allowed to mate randomly for 5 generations. We chose 

this approach, as a compromise to minimize confounding laboratory adaptation while 
ii2 still incorporating genetic variation present in natural populations. With this approach 

linkage disequilibrium among variants will likely be more extensive than in wild-caught 
H4 flies. Post-intercross, we maintained the population at large size on cornmeal molasses 

media with live yeast in four 200ml culture bottles at 24°C and 60% humidity. 
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lie Predation Environment 

We used first instar nymphs of the Chinese mantid (Tenodera aridifolia sinensis) as 
us predators. We collected mantid egg cases locally from old fields in southern Michigan 

and supplemented these from garden suppliers (Nature's Control Medford, Oregon) when 
120 necessary. We hatched and maintained egg cases at 24°C and 60% humidity. Approxi- 
mately 100-400 mantids emerged from each egg case and were used as predators for the 
122 duration of the first nymphal instar. After hatching, we housed mantids at 18°C and 60% 

humidity in arenas consisting of a 710mL plastic cup with a mesh covered window for air 
124 flow. We placed a tissue at the bottom of each cup to trap moisture when watering to help 

maintain humidity. We also added a green plastic aquarium plant to provide substrate 
126 for mantids to perch upon. Unless otherwise specified we used five mantids per arena. 

All episodes of predation occurred at 18.5°C and 60% humidity, and were initiated 
128 between 12-3 pm. We fasted mantids for 24 hours before each episode of selection to 

increase predation rates. Arenas were cleaned with 70% ethanol and water before use. 
130 25 flies were introduced into each arena via a funnel, after which arenas were returned to 

the incubator. After 24 hours, all predation arenas were moved to a 4°C refrigerator to 
132 knock down flies and mantids to aid collecting. We then removed the mantids from each 

container, and surviving flies were censused. 

134 Testing the role of flight in the escape response using a vestigial 1 
mutant population 

136 To test whether flight played a role in the escape response, we tested whether wing loss 
would negatively impact survival under risk of predation. We introgressed a mutation in 
138 the vestigial (vg 1 ) gene into the base population by repeated backcrossing for 10 gener- 
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ations.The vg 1 mutation causes a nearly complete loss of the wing blade and associated 
140 musculature (Sudarsan et al., 2001). We competed individuals from the vg 1 mutant pop- 
ulation with their wild-type conspecifics by placing 13 mutant and 13 wild-type flies in 
U2 each of 16 arenas (8 arenas each for male and female flies). The survivors for both vg 1 
and their wild-type conspecifics were counted after 24 hours with the predators. 

144 Assaying Phenotypic Selection: Base population 

To assess how naturally segregating variation for wing size and shape might be associated 
146 with survival during predation events, flies from the base population were exposed to 

predation. Predation on males and females was assayed in separate arenas so that we 
148 could examine independent effects of sex. We placed 20 flies into each arena (9 arenas 

each for females and males). Four days later we set up a second block of arenas (10 arenas 
150 for females, 8 for males). After predation, we collected all surviving flies and all wings 

from consumed flies from the bottom of the arenas. We also collected 100 individuals 
152 of each sex that were not exposed to predators. All bodies and wings were preserved in 

ethanol for dissection and measurement. 

154 Experimental Evolution 

We randomly selected five hundred flies from the base population and used these as par- 
156 ents to generate the four populations for experimental evolution. We randomly assigned 

offspring of these 500 individuals to each of the treatments, with blocks of offspring for 
158 the different replicates. The predator free populations control for selection and adaption 

independent of the predators (i.e. alterations in competitive environment). Selection 
i6o was administered in two replicate sets each consisting of one predation and one control 
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population, hereafter referred to as PredRl, ConRl, PredR2, and ConR2. We offset the 
162 generational cycle of replicate 2 from replicate 1 by 2-5 days for logistical reasons, but 

the populations were otherwise treated identically. Each population was reared in four 
164 bottles with approximate densities of 100-500 eggs per bottle each generation. We did not 

explicitly control for density, but restricted egg laying time to 2-6 hours to avoid larval 
lee overcrowding. Bottles were reared at 24° C and 60% humidity until eclosion of adults. 

Three days after eclosion of the first flies, progeny from each population were lightly 
168 anesthetized using C0 2 , placed randomly into vials, and maintained at 18.5°C and 60% 

relative humidity. The following day, flies from a given treatment were mixed under 
no anesthesia, to minimize inadvertent selection on developmental time. Flies (25/vial) were 

transferred to fresh vials, corresponding to the number of arenas used for predation in 
172 that generation. Each generation we varied the total number of arenas depending on 

the voracity of the current batch of mantids in order to ensure that the total number 
174 of survivors was large enough to limit the effects of drift (between 150-400 surviving 

individuals/generation). Flies from the control vials were similarly mixed under anesthesia 
176 after which we placed 50 flies into each of 8 vials. Remaining predation and control flies 

were set aside as backups. All flies were then returned to the incubator at 18.5°C and 60% 
178 for at least 24 hours prior to the episode of predation. Because egg cases were seasonally 

available, we occasionally used second instar mantids to maintain experimental evolution. 
i8o However, second instar mantids were not used for any experimental trials. Control arenas 

were identical to predation arenas, only lacking mantids. 
182 After selection, we collected all survivors from the predation arenas. To maintain 

similar population sizes we selected individuals at random from the control populations 
184 matching the number of male and female survivors from the respective predation pop- 
ulation. Individuals from each treatment were transferred into separate 30x30x30 cm 
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186 polyester mesh cages, and allowed to recover for 30-45 minutes before fresh bottles of 
food media were placed into each cage. After allowing sufficient time for egg laying, the 

188 bottles were then removed from the cages and reared at 24° C and 60% humidity After 
breeding, remaining adult flies were stored in ethanol at -20°C. 

wo Assaying Phenotypic Selection: Evolved populations 

To examine how the fitness function changed as a result of experimental evolution, we 
192 repeated phenotypic selection (as described above) during generations 31 and 32 of the 

experiment. Given the large size of this experiment, it was performed in four blocks, 
194 with two blocks for each generation. At generation 31 of experimental evolution, we set 

up 14 arenas each of PredRl females & males and ConRl females & males. Five days 
we later, we set up 14, 14, 8, and 9 arenas for PredR2 females & males and ConR2 females & 

males respectively. At generation 32, we set up 14 arenas each of PredRl females & males 
198 and ConRl females & males . Five days later we set up 13, 13, 14, and 14 arenas for 

PredR2 females & males and ConR2 females & males respectively. As before, we collected 
200 all surviving flies and unconsumed wings and stored them in ethanol. Overlapping egg 

cases were used for this experiment, and egg case of origin was used as a covariate in the 
202 model (see below). We distributed mantids so as not to confound predation effects across 

replicates and treatments. 

204 Wing Measurement &; Statistical Analysis 

Wings were dissected and mounted on slides in 70% glycerol. When available, both 
206 wings from an individual were mounted. Wings were also dissected from 25 flies that 
were stored from the initial generation of experimental evolution, and from every 10 



Downloaded from http://biorxiv.org/on September 18, 2014 



11 

208 generations following up to generation 30 to estimate the trajectory of size and shape 
change. Wings were imaged at 40X magnification on an Olympus DP30BW camera 

210 mounted on a Olympus BX51 microscope using 'DP controller' V3.1.1 software. All 
images were saved in greyscale as tiff files. 

212 To capture landmark and semi-landmark data we followed a modified protocol (Pitch- 
ers et al., 2013) for the use of the wingmachine software (Houle et al., 2003). We used 

2u the program tpsDig2 (Rohlf 2011) to manually record the coordinates of two starting 
landmarks, and used wingmachine to fit nine B-splines to the veins and margins of the 

216 wings in the images. We extracted 14 landmark and 34 semi-landmark positions, and 
performed Procrustes superimposition (Zelditch et al., 2012). After superimposition, the 

2i8 positions of semi-landmarks were allowed to slide along each segment of the wing mar- 
gin/veins, minimizing Procrustes distance, using CPR V0.2 (Marquez 2011). The data 

220 were checked for visual outliers at multiple stages; and putative outlier images were reex- 
amined and splines re-fit if necessary. The (semi-) landmark configurations for all wings 

222 measured for this study were superimposed together, resulting in a common shape space. 
Centroid size was used as measure for size (Zelditch et al., 2012). For flies with both 

224 wings collected, we calculated the mean shape and centroid size per individual. 

Model Selection 

226 For the univariate analyses, we evaluated model fits using Akaike's Information Criteria 
(AIC) and Bayesian information criteria (BIC). AIC has been shown to often 'prefer' 

228 more complex models than there is support for, particularly when sample sizes are large 
(Grueber et al., 2011), we used model weights from BIC throughout for consistency to 

230 perform model averaging when appropriate. Unless otherwise noted, all further analyses 
were conducted in R (V2.15.1) (R Core Team, 2012). Scripts for the custom functions 
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232 described below are available with the data from DRYAD. doi.XXXX. 

234 Analysis of Survival 

For the vg 1 mutant and wild-type competition assays, we fit the model: 

WT prop ~ N{fi + /3 sex ,a 2 ) (1) 

236 where WT prop was the proportion of wild-type survivors in each arena, and (3 sex was the 

model coefficient for sex. 
238 For the base and evolved populations, we measured survival ability as the total number 

of surviving flies in each container. For the base population, we fit the model: 

Surv ~ N(ii + p sex ,a 2 ) (2) 

240 along with a set of expanded and restricted models (Supplementary table 1) where Surv 
was the number of surviving flies in each arena and (3 sex was the model coefficient for 

242 sex. Model averaging produced coefficient estimates indistinguishable from the model 
with best support so this model was used for further inference. 

244 For the evolved populations we fit the model: 

Surv ~ N(fl + P S R + Pgen + Peggcs, ^) (3) 

along with a set of expanded and restricted models (supplementary table 2) where Surv 
246 was the number of surviving flies in each arena and /3sr, (3 ge n, and f3 eggcs were the model 

coefficients for selection regime, generation of selection when the assays were performed, 
248 and the egg case of origin for the mantids in each arena. The coefficient estimates produced 
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from model averaging models 1 and 2, which accounted for 95% of the overall support, 
250 were indistinguishable from the model with best support so for simplicity we used it for 

further inference. For the above models, we confirmed the effects using generalized linear 
252 models (poisson with log link), or a logistic regression with similar results. 

Analysis of Phenotypic Selection on Size 

254 We used the Lande and Arnold (1983) approach to examine selection acting on wing size 
in the base and evolved populations. As recommended by Janzen and Stern (1998), we 

256 used logistic regression on survival for statistical inference and general linear model on 
relative fitness to estimate coefficients. Relative fitness for each individual was calculated 

258 by scaling survival (0 for dead and 1 for survived) by the total proportion of survivors in 
each experiment. To measure linear selection in the base population we fit the model: 

pr(W) ~ 6z(p = logir\0o + /W) (4) 

260 along with a set of expanded and restricted models (supplementary table 3). W was ab- 
solute fitness (survival) and fi S i ze was the model coefficient for standardized wing centroid 

262 size. The coefficients produced by model averaging were indistinguishable from the model 
with best support so it was used for further inference. It should be noted that we are 

264 estimating the linear S, and non-linear C selection differentials (Brodie III et al., 1995). 
The /3's in the equations are used to represent estimated model parameters, and do not 

266 represent selection gradients. 

To measure linear selection in the evolved populations we fit the model: 
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pr(W) ~ (p = 

logif {fio + P size + PsR + Prep + P sex Pgen PsizexSR Prepxgen Psexxgen + PsRxrep)) 

(5) 

268 along with a set of expanded and restricted models (supplementary table 4) where W 
was fitness and P S i ze , Psr, Prep, Psex, and P gen were model coefficients for standardized 

270 wing centroid size, selection regime, replicate, sex, and the generation of experimental 
evolution respectively. We fit separate models as above including the quadratic effect of 

272 size to estimate non-linear selection. Estimates for non-linear selection on size were near 
zero, non-significant, and did not improve model fits in either the base or the evolved 

274 populations. 

Non-parametric estimation of the form of the fitness functions substantially aids vi- 
276 sualization and interpretation of fitness functions (Schluter, 1988). We therefore used 

generalized additive models from the the MGCV package VI. 7. 22 (Wood, 2004) to fit cu- 
278 bic splines to subsets of the data from each experiment corresponding to the relevant 

significant effects estimated by the logistic regression analyses. Optimal smoothing pa- 
280 rameters were estimated using REML. 

Variances in Size & Shape 

282 One additional approach to investigating natural selection is to examine the changes in 
phenotypic variance before and after the selective event (Endler, 1986). Under either 

284 directional or stabilizing selection, a reduction in variation would be predicted. Under 
disruptive selection however, we would predict an increase in variation. Analyzing dif- 
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286 ferences in variance between dead and surviving flies, as well as between predation and 
control populations may therefore provide additional information on the type of selection 

288 occurring in these populations. We used Levene's test to assess changes in variance for 
wing size, using deviations from the median rather than the mean since this approach is 

290 more robust to departures from normality. For the base population, we modeled the main 
effects of sex and size because our previous analyses lacked support for an interaction 

292 between sex and the form of selection. We fit the model: 

Ld ~ N(ji + P sex + Pw, o 2 ) (6) 

where Ld was the Levene's deviates for each individual, (3 sex was the model coefficient for 
294 sex and Pw is the model coefficient for absolute fitness. Though our previous analyses do 

not suggest that selection acting in the evolved populations differed between replicates, 
296 differences in size between PredRl and PredR2 suggest that its inclusion is appropriate. 

We fit the model: 

Ld ~ N(jl + (3 S R + Prep + Psex + Pw + PsRxrep, & 2 ) (7) 

298 where Ld was the Levene's deviates for each individual and fisR, Prep, ft sex, and Pw were 
the model coefficients for selection regime, replicate, sex, and absolute fitness, respectively. 

300 Confidence intervals for all estimates were generated by non-parametric bootstraps, in 
order to avoid issues with non- normality of residuals. We also calculated the coefficient of 

302 variation for each of the groups modeled above. Because we know that there is substantial 
sexual size dimorphism as well as size differences among the base and evolved populations, 

304 this approach should provide more intuitive visualization. 

To compare levels of variation in shape we took a somewhat simpler approach. We 
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306 expressed the variability of each group as the trace of its covariance matrix for shape. 

We then bootstrapped the data to generate samples of each covariance matrix in order 
308 calculate confidence intervals on the estimated matrix trace. Non-overlapping (95%) 

confidence intervals intervals were then used to infer statistical support for differences in 
310 variance among groups. 

Multivariate Analysis of Shape 

312 In our initial analyses, we found that the modeled effects of allometry and sexual dimor- 
phism were extremely consistent between treatments and over time (i.e. the vectors of 

314 model coefficients for sex and for size were very tightly correlated; see below). In order to 
facilitate the interpretation of the modeled coefficients of selection and generation num- 

3i6 ber, we therefore sought to exclude these effects from our analyses. With data from all 
wings pooled, we fit the model: 

S ~ N(p + {3 sex + {3 size , S) (8) 

3i8 where S is the matrix of Procrustes coordinates and f3 sex and (3 size are the vectors 
of model coefficients for sex and for wing centroid size respectively. E is the "error" 

320 covariance matrix. We retained the residuals from this model as our shape variables. 

Configurations of Procrustes coordinates by definition include dimensions without vari- 

322 ance. The Procrustes superimposition results in a deficiency of 4 ranks (1 each for removed 
size and rotation information, and 2 for position), and each semi- landmark may contribute 

324 as little as 1 added dimension (Zelditch et al., 2012). In order that the shape data would 
not be rank deficient, we extracted principal components from the (96-dimensional) resid- 

326 uals, and retained the first 96 — (4 + 34) = 58 principal components, comprising > 99.9% 
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of the shape variance in the full set of residuals. Shape PC's used in all the analyses below 
328 are thus of full rank, and are expressed in a common sub-space. 



330 Modelling Shape Change 

Separately within each of the four evolved populations, we estimated the direction of 
332 observed evolutionary change as the vector of model coefficients from the multivariate 
linear model: 

S p ~iV(^ + /3 5en ,S) (9) 

334 where S p is the matrix of principal component scores for shape in a given population 
and f3 gen is vector of model coefficients for time, expressed as the number of generations 

336 removed from to the base population. Once we had estimated these vectors of parameters, 
we compared their directions by calculating vector correlations as: 



la ■ b| 



(10) 



||a|| x ||b|| 

338 where |a- b| is the absolute value of the dot (scalar) product between vectors a and b, 

while ||a||, and ||b|| are the magnitudes (L 2 , or Euclidean norm), for each vector. The 
340 absolute value for the dot product was used to avoid any numeric issues with arbitrary 

sign changes that can occur computationally (during the bootstrapping procedure, see 
342 below). Thus r = 0 represents no similarity between the vectors while r = 1 means the 

two vectors point in an identical orientation (but possibly opposite in direction). Given 
344 that r is a multivariate extension of the Pearson correlation co-efficient p, we consider this 

a more intuitive measure than the vector angle (8 = arccos (r) in radians) which has been 
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346 used elsewhere. Confidence intervals were computed using non-parametric random pairs 

bootstrapping, from 10,000 bootstrap iterations. This approach was used both to compare 
348 the direction of S as measured in all five populations, and to compare the directions of 

observed shape change among the evolved populations. 
350 To illustrate the magnitude of change in wing shape during experimental evolution, 

we calculated a shape score (Drake and Klingenberg, 2008). Briefly, we projected the 
352 shape data onto a line in the direction defined by the vector of model coefficients for the 

generation term (/3 gen ) from model (3): 

shapescore = Y/3 T (/3/3 T )" a5 (11) 

354 The shape score provides a univariate measure of shape change that can be plotted 
against generation number to visually assess the magnitude and linearity of the rela- 

356 tionship (Drake and Klingenberg, 2008). We used custom R functions to calculate vector 
correlations and shape scores. 

358 Selection on Shape 

Within each population, we estimated the vector of linear shape differentials (S). Tra- 
360 ditionally this would be calculated as vector of differences between the mean phenotype 

of survivors and the mean phenotype of those individuals that were preyed upon. Here 
362 we estimated S using a 2-block partial least squares (PLS) approach (Rohlf and Corti, 

2000; Klingenberg and Zaklan, 2000; Mitteroecker and Bookstein, 2011; Klingenberg and 
364 Monteiro, 2005; Gomez et al., 2008, 2006) with the matrix of the 58 shape PC's forming 

one block, and the vector of survival data (0 or 1 for dead or survived) as the second block. 
366 We note that in this case this estimate of S is proportional to the differences between the 
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mean shape configurations for the dead and survivors. 

368 It is important to note that wing shape itself is the trait, and not individual land- 
marks/PC's. After Procrustes superimposition, individual landmarks and semi-landmarks 

370 cannot be meaningfully interpreted independent of the whole shape configuration and the 
superimposition can generate correlation between landmarks that is confounded with bi- 

372 ological correlations (Zelditch et al., 2012). Thus, interpreting the selection gradients, 
f3, from a multiple regression for shape (sensu Lande and Arnold, 1983) for "individual" 

374 shape variables is biologically meaningless (Albert et al., 2008). In addition, selection 
gradients can be difficult to visualize for shape (Klingenberg and Monteiro, 2005) but see 

376 (Mitteroecker and Bookstein, 2011), in particular because estimating the inverse of the 
phenotypic covariance matrix, P _1 , can be problematic. We observed that, upon resam- 

378 pling, lack of stability in P _1 caused computationally difficulties. One alternative is to 
retain only the first few PC's and analyze them as if they were independent traits (Gomez 

380 et al., 2006, 2008; Kuchta and Svensson, 2014). This is still sub-optimal, however, since 
substantial variation and selection may be missed and the biological interpretation of any 

382 selection that is detected is difficult. While this is an important and outstanding issue, we 
elected to use selection differentials for the shape analyses because they retain biological 

384 meaning and the focus of the study is on the predictability of selection not its specific 
form or estimation. However, this does mean that the results need to be interpreted as a 

386 combination of both direct and indirect selection on shape. 

We estimated total selection on wing shape as the magnitude of the vector of the 

388 selection differentials, ||S||, and used sampling with replacement of the data to gener- 
ate non-parametric bootstrap confidence intervals on these estimates. Additionally, we 

390 permuted survivorship relative to the measures of shape to assess the null hypothesis 
that wing shape does not contribute to variation for survivorship. We also compared 
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392 the directions of the S vectors using vector correlations as described above. Finally, we 
wanted to assess the degree to which the experimental evolution populations had evolved 

394 in the direction 'predicted' by selection as measured in the base population. To do this 
we calculated the vector correlations between the S vector measured in the precursor 

396 population and the vector of model coefficients for generation (/3 gen from model (8)) as 
modeled separately for each population. 

398 Results 

Evidence that flight aids in the escape response 

400 To test whether flight performance and wing form were potential targets of selection driven 

by the mantid predators, we introgressed a mutation in the vestigial (vg) gene into our base 
402 outbred population that nearly completely ablates the wing blade and associated flight 

muscles. We competed vg 1 (functionally wingless) flies and their wild-type conspecifics 
404 with the predators. As predicted, the vg 1 individuals were disproportionately the targets 

of predation. The survivors for both sexes consisted of approximately 60% wild-type and 
406 40% mutant individuals (figure SI), consistent with a role for flight and possibly wing 

morphology in the escape response of Drosophila. 

408 Predator driven selection on natural variation for wing form 

We next asked how natural variation for wing form was associated with survivorship 
4io by exposing flies from the base population to the mantids. We observed evidence for 

significant negative directional selection on wing size (Figure IB, Table SI, S = -0.29 
412 ±0.22, p ~ 0.01) with little evidence for non-linear selection (C = -0.05 ±0.22, p ~ 

0.16). Visualization by fitting cubic splines to the survival data (Schluter 1988) was 
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4u consistent with the estimates of directional selection (figure IB). Despite sex specific 

differences in survivorship (6.7 ±0.75 & 3.8 ±1.07 survivors per arena for females and 
4i6 males respectively), evidence was weak for an interaction between selection on size and 

sex (Table SI). It is currently unclear whether the target is wing size per se or whether 
4i8 it is due to a correlation between wing size and overall body size. 

Additionally, shape has been shown to be correlated with escape ability in other or- 
420 ganisms (Langerhans et al., 2004; Dayton et al., 2005; Langerhans and Makowicz, 2009; 

Svensson and Friberg, 2007). Using a 58 dimensional representation of wing shape (figure 
422 S2), we used partial least squares (PLS) to estimate the vector of selection differentials 

for shape (S) in the base population. We visualized this S vector for shape comparison 
424 (figure 1C). We see selection for a change in aspect ratio in which relatively longer and 

narrower wings are favored. 

426 What does wing form look like after experimental evolution? 

We asked whether experimental evolution under risk of mantid predation would result 
428 in phenotypic changes in wing form consistent with our initial estimates of selection. 

Previous work has demonstrated that both wing size and shape have moderate to high 
430 heritabilities, and most aspects of wing shape are readily altered by artificial selection 

(Mezey and Houle, 2005; Weber, 19906). To assess this we subjected two populations 
432 of flies to episodic selection by mantids, each paired with control populations evolved 

without predators. The control populations allow us to assess the effects of selection due 
434 to the experimental evolution procedure unrelated to the predators that occurred during 

the experimental evolution process. 
436 As expected, survival in both predator populations increased, compared to the controls 

(10.86 ±0.78 & 13.98 ±0.79 survivors per arena for control and predation populations re- 
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438 spectively; figure 2A). This represents a ~30% increase in survivorship relative to the 

control populations. We did not observe differential survival between sexes in this exper- 
440 iment for either selection regime. 

We measured individuals stored during the experimental evolutionary process every 
442 ten generations, from the base to generation 50, to track changes in wing form. Wing 

size of all populations increased ~3.7% over the 50 generations of experimental evolu- 
444 tion (0.0035 mm ±0.0004 mm per generation). This change was most likely a result of 

selection due to non-predatory aspects of the experimental evolutionary procedure be- 
446 cause all populations increased at similar rates and maintained the same overall sexual 

dimorphism. However, environmental variation in these samples collected directly from 
448 the experimental evolution regime is relatively large (Figure S3). 

To more carefully estimate size differences among the evolved populations, we mea- 
450 sured wing size in the overall population by using the wings from the dead and surviving 

flies from the phenotypic selection experiments, as all flies were reared under density con- 
452 trolled conditions. Thus environmental and genetic effects were not confounded. Com- 
parison of the number of surviving and dead flies from this assay to the number of wings 
454 recovered suggests that nearly all wings from dead individuals were recovered, and should 

provide reasonable estimates. 
456 Under these conditions, the relevant contrast is the difference between the control 

populations and predation populations. We found that the two control populations had 
458 similar wing sizes, yet the two populations evolved under risk of predation diverged in 

size (figure 2B) even though all populations showed a general size increase relative to the 
460 ancestral population. Surprisingly only PredR2 has evolved in the predicted direction, 

with wings ~1% smaller than controls, while PredRl evolved wings that are ~2% larger. 
462 In terms of shape, all four populations have evolved from the base population, though 
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not to an equal extent. We visualized the evolutionary trajectories of the four populations 
464 by plotting the shape score for generational effects (equation 10) (figure 3A & B). In all 

four cases the evolutionary trajectories were best described by a simple linear model. 
466 Whereas the two control populations have changed in a very similar fashion (figure 3A), 

the two predation populations are clearly divergent, with wing shape in PredR2 evolving 
468 significantly more rapidly than in PredRl (figure 3B). 

Over the course of experimental evolution the wings of all populations have changed 
470 aspect ratio: their length increasing slightly as their depth decreases. This change is most 

pronounced in PredR2 (figure 3C). Other than the differences in aspect ratio, PredRl 
472 and PredR2 differ most noticeably in the response of the cross-veins and the distal end 

of L5. PredRl demonstrates a proximal shift in the posterior cross-vein and an anterior 
474 shift in the attachment of L5 to the margin; by contrast in PredR2 there was no change 

in L5 and an anterior shift in the anterior cross- vein. 

476 What do the fitness functions look like after experimental evolu- 
tion 

478 After 30 generations of experimental evolution, we again exposed flies that evolved with 

(and without) predators to a bout of predation. We observed negative directional selection 
480 for size in the control populations (S = -0.16 ±0.07, p ~ 0.0001; figure 4A), consistent with 

our findings from the ancestral base population. We also observed negative directional 
482 selection in the predation populations, but of diminished magnitude relative to the controls 

(S = -0.06 ±0.05, p ~ 0.0005; figure 4B). Both the control and the predation populations 
484 showed extremely weak quadratic selection on size (C = 0.0001 ±0.03, p ~ 0.73). We 

might reasonably expect that the reduction in the magnitude of directional selection in 
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486 PredR2 was a result of evolutionary change in wing size in response to selection. However, 
this provides no explanation for the increase in size in PredRl. 

488 We assessed selection on wing shape in the evolved populations as S; the vector of 
selection differentials between captured flies and survivors and compared the magnitudes 

490 of total selection on shape from the differential, ||S||. For the estimates of ||S||, we also 
generated distributions under the null expectation (of no association between wing shape 

492 and survival) using permutations of the data. In addition we also calculated the vector 
correlations between differentials in order to quantify their degree of alignment. For 

494 both approaches we computed confidence intervals on our estimates by applying a non- 
parametric bootstrap approach. 

496 As can be seen from figure 4C, there was evidence for a significant association be- 
tween shape and survival in the presence of the mantid predators for all populations: 

498 the estimates of ||S|| exceeded the 95% threshold permuted under the null hypothesis. 
Also notable are the much smaller ||S|| estimates of the predation populations compared 

500 to that in the base population. This evidence is consistent with a relative reduction 
in the magnitude of selection experienced by the predation populations after 30 genera- 

502 tions. Interestingly, there is some evidence for difference in ||S|| between the two control 
populations, however both still exceed the predation treatment regimes. 

504 The directions of the S vectors in the predation populations are also quite divergent, 
however, and their vector correlations with the base differ; 0.39 (0.07 - 0.59) & 0.02 (~0 - 

506 0.34) though our estimate are not very precise (vector correlation followed by bootstrapped 
95% confidence intervals, base S vs. PredRl S & PredR2 S respectively). By contrast 

508 the vector correlations between the controls and the base population were consistent with 
one another at 0.39 (0.01 - 0.51) & 0.29 (0.04 - 0.48) respectively (see figure S4 for all 

510 pairwise comparisons). Moreover, the directions of the S vectors measured in the four 
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evolved populations are not closely aligned (figure S5). 

512 Has wing form evolved in the direction predicted by selection on 
the base population? 

514 While there has been evolution of shape in all four populations, we wanted to assess how 

much of the observed change is in the predicted direction (based on S in the ancestral 
5i6 base population). We calculated the vector correlations between the generation shape 

change vectors from each evolved population and the S vector from the base population 
5i8 and observed that the evolutionary responses of the predation populations were more 

aligned with the predicted vector compared with the control populations (figure 5). It 
520 is notable that none of the populations are particularly highly aligned with the initial 

predicted vector. 

522 Given that both predation populations have experienced a similar reduction in the 
magnitude of selection (as represented by ||S||: figure 4C), and a similar amount of evo- 

524 lutionary change in this direction (figure 5), it appears that they experienced similar 
changes in the selective function for shape, despite their divergent evolutionary response 

526 for size (figure 3C). This suggests that the two predation populations are evolving dif- 
ferent avoidance strategies in response to the predation pressure imposed by the mantids 

528 - likely involving traits other than wing morphology. In the case of PredR2 there is 
evidence that the reduction in the intensity of selection is associated with evolution in 

530 the predicted direction, but it seems likely that there may be other adaptions occurring 
in PredRl 
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532 Changes in variance in both size and shape 

Changes in variation can also provide information about the form of selection experi- 
534 enced by a population. Evidence for differences in variance for size between survival 

classes in the base population was weak (Levene's deviates = 0.28 ±0.04, 0.27 ±0.03, 0.24 
536 ±0.04 for unselected, dead, and surviving individuals respectively, p ~ 0.25). Males had 

lower variance for all survial classes ( -0.07 ±0.03, p < 0.001). For the evolved popula- 
538 tions, differences in size variance were only found among populations, with ConRl and 

PredR2 having equal variance (Levene's deviates = 0.11 ±0.01, 0.11 ±0.01 respectively), 
540 ConR2 having higher variance (Levene's deviates = 0.15 ±0.01), and PredRl having lower 

variance (Levene's deviates = 0.10 ±0.01). Surviving individuals trended towards lower 
542 variance, but the difference was not significant (-0.003 ±0.006, p ~ 0.37). Males again 

had lower size variance, but a much lower magnitude of difference (-0.015 ±0.006, p < 
544 0.005). 

Estimates of variance for shape in each population (the trace of the covariance matrix) 
546 show a dramatic reduction in variance in surviving flies and lower overall shape variance 
in the populations that evolved under predation risk (figure 6). Not only is the variance 
548 lower in the predation populations as compared to the controls for shape, but the sur- 
viving populations have much lower variation when compared to the populations that 
550 were captured and eaten by the mantids suggesting that selection has already reduced 
variation in the evolved populations and continues to do so. 

552 Discussion 

For several decades, phenotypic selection analysis has been used to attempt to identify 
554 the primary targets of selection within natural populations under the assumption that the 
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presence of selection on specific traits would provide information about about how those 
556 traits evolved and what future changes could be expected. The striking levels of conver- 
gence and parallellism in several well known study systems suggests that this assumption 
558 may be valid (O'Steen et al., 2002; Langerhans et a!., 2004; Dayton et al., 2005; Langer- 

hans and Makowicz, 2009; Auld and Relyea, 2011; DeWitt et al., 2000, 1999). However 
560 estimates of phenotypic selection analyses have not always been able to predict long-term 

evolution. It remains unclear whether this contradiction is a result of publication bias 
562 or indicative of larger issues. In this study we measured phenotypic selection on a naive 

population in response to a novel predator. We then re-measured the strength and direc- 
564 tion of natural selection after populations were allowed to evolve under natural selection 

with the predator to determine whether our results would match the patterns of parallel 
see response cited above. We found that the populations evolved divergent morphology for 

size but relatively consistent shapes. What do these results tell us about the form of 
568 natural selection, and what are the implications for its use in evolutionary prediction? 

To use the breeder's equation to predict an evolutionary response, we require not only 
570 a vector of directional selection, but heritable variation along the same axis as selection 

(Hine et al., 2011; Walsh and Blows, 2009). The direction of this genetic variation is 
572 determined by the size and structure of the genetic covariance matrix G. A number 

of studies have demonstrated that populations tend to evolve along genetic lines of least 
574 resistance (Schluter, 1996; Mcguigan et al., 2005), not necessarily the direction of strongest 

selection. 

576 Previous work has demonstrated that there is considerable segregating genetic vari- 
ation in most populations for wing shape. In particular the effective dimensionality of 

578 G for wing shape is quite high, and close to the number of measured traits (Mezey and 
Houle, 2005). Given the size and scope of the endeavor, we did not attempt to estimate 
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580 G for the base population we used. However, genetic variation among the progenitor 

strains used to generate the population shows a high effective dimensionality (data not 
582 shown), consistent with previous results from other populations. It is possible that genetic 

variation in the direction of selection imposed by mantid predation may be minimal, and 
584 that the genetic line of least resistance is not perfectly aligned with this direction. Thus 

at least some of the common changes in wing form may be the result of a combination 
586 of lab domestication and evolution along the genetic lines of least resistance. Despite 

this, we see clear evidence for more shape change in the predation regimes consistent 
588 with the initial vector of selection. Given the high dimensionality (58) of shape, this is a 

pronounced effect, demonstrating that even with potentially countervailing selective and 
590 genetic forces, selection is still altering shape as predicted. 

In addition to the need for available genetic (co)variation, there are several factors 
592 that can influence the evolutionary response to directional selection including: indirect 

selection due to correlated traits, stabilizing selection, fluctuating selection, fitness trade- 
594 offs, and unmeasured (but strongly selected) traits (Kingsolver and Diamond, 2011). The 

selection differentials reported include direct and indirect selection, so even though we 
596 cannot estimate the separate contributions of each, we saw significant total selection on 

both wing size and shape. Additionally, neither the base nor the evolved populations 
598 showed evidence for stabilizing selection that might have reduced the strength of direc- 
tional selection. Though we expected the base population to be far from any fitness peak, 
6oo it is curious that the predation populations did not show stronger non-linear selection 

given the reduction in the strength of the linear component. 
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602 How similar is the form of selection after 30 generations of ex- 
perimental evolution 

604 The selective pressure imposed by the mantids on size remained relatively stable, con- 
sistent with the conclusions reached by Morrissey and Hadfield (2012). Our measure of 

606 S in the ancestral base suggested strong directional selection for smaller wings (~-0.29). 
After nearly two years and 30 generations of evolution, the estimate of S was remarkably 

608 similar in the control populations (-0.16, figure 4). This slight reduction in the magnitude 
is, perhaps, unsurprising since the difference in sample size — nearly an order of magnitude 

6io greater for the evolved populations — allowed for more precise estimation. The estimate 
is also in line with the median reported by Kingsolver and Diamond (2011) for size traits 

612 (|0.14|) and viability (|0.08|). As a result, there is little evidence to suggest that temporal 
variation in selection from generation to generation resulted in the divergence in size in 

6u the predation populations particularly because both populations still exhibit measurable 
directional selection. 

6i6 Selection on shape 

For shape, the picture is less clear. Because of the high dimensionality of shape, not only 

6i8 is estimation much more difficult, but there is a much larger available phenotype space. 
Perhaps unsurprisingly then, the vector correlations between selection in the base popu- 

620 lation and the control populations are reasonably low (~0.35) suggesting that fluctuating 
directionality of selection may have limited the evolutionary response. The degree to 

622 which there was true variation in direction of selection for shape, as compared to estima- 
tion issues (even with our large sample sizes) remains unclear. Indeed, this is one of the 

624 major reasons we used S instead of (3, as estimating P -1 proved to be computationally 
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difficult, and caused problems during resampling. Despite this, both PredRl and PredR2 
626 show considerable overlap between the vector of shape change during evolution and the 

direction of selection predicted in the base population (r ~0.5, figure 5). This suggests 
628 that even though the form of selection for shape is apparently less stable than for size, it 

has not resulted in substantial divergence between populations. 
630 It is worth considering what is lost by using the selection differential S instead of the 

gradient, (3 = P _1 S. For most phenotypic selection studies the main difference relates 
632 to disentangling direct and indirect selection on traits (pre-multipling by P _1 removes 

the phenotypic covariation). Shape data is unique, in that the different variables are not 
634 independent traits. Instead the whole configuration (as represented by a vector for each 

individual) is a geometric representation of the shape "trait". Pre-multiplication by P 1 
636 has the potential to change the observed orientation of the vector of the selection differ- 
ential S, however it also causes difficulties with interpretation of the resulting selection 
638 gradients, and so the preferred method is to visualize the selection differentials (Klingen- 

berg and Monteiro, 2005) as we have done here (but see (Mitteroecker and Bookstein, 
640 2011) for an alternative perspective). Other groups have instead utilized a small number 

of principal components of the shape data in a standard Lande-Arnold selection gradi- 
642 ent analysis (Gomez et a!., 2006; Kuchta and Svensson, 2014). However, this utilizes a 

fraction of the variation in shape, with no guarantee that it represents the components 
644 of variation under selection. Thus a full multivariate approach is needed (Klingenberg, 

2010) though we currently lack an accepted standard sensu Lande-Arnold . We suggest 
646 that continued effort and discussion into estimating and visualizing selection on shape, as 

well as determining the appropriate "dimensionality" of such effects is warranted. 
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648 Possible causes of divergence and parallel evolution 

Unknown fitness trade-offs and lab adaptation may have played a role in the divergence 
650 between PredRl and PredR2. During the course of evolution, all four evolved populations 

showed a net increase in wing size of ~3.2 % and a lengthening and broadening of the wing 
652 blade in direct contrast to the smaller, longer, and narrower wings favored by selection in 

the base population (figure 1). These changes were remarkably consistent among the four 
654 populations and are likely a result of selection due to shared aspects of the experimental 

evolutionary process independent of the predators. However, though the directionality of 
ese the shared evolved response and selection measured in the base population suggests that 

evolution may be slowed in the predation populations, this gives little indication as to 
ess why PredRl diverged from the predicted size trajectory. This could be caused by drift 

between the replicates over the 30 generations of experimental evolution. However, this 
660 explanation is unsatisfactory as the evidence of drift is missing in the control populations 

wherein it would be expected to dominate. Though the utmost care was taken to control 
662 variation between the replicates, differences in the health and voracity of the mantids, 

as well as in some other environmental factors was unavoidable, possibly contributing to 
664 this effect. 

Where do these results leave us? We possess robust theory for measuring selection, 
eee and for predicting evolutionary responses into the near future (Lande and Arnold, 1983). 

However, we are often left to assume that populations will evolve phenotypes in the distant 
668 future consistent with these estimates. Though a number of other systems have examined 

the evolutionary consequences of manipulating predation regimes long term, notably the 
670 work of David Reznick and colleagues (Reznick et al., 1990, 1997; Reznick and Ghalambor, 

2005), few studies have investigated how well evolutionary responses coincide with specific 
672 measures of selection. The Drosophila-mantid system described here allows us to maintain 
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specific selective pressure in a relatively homogeneous environment on a population with 
674 a known history. This allows us to not only impose specific selection pressures, but to 

remeasure selection itself during the evolutionary process. 
676 It is likely that unmeasured anti-predator behavioral traits played an important role. 

Unmeasured traits that may be under selection (and genetically covary with measured 
678 traits) can profoundly influence the biological inferences we make about natural selection, 

and evolutionary response. While many studies of phenotypic selection attempt to ex- 
680 amine multiple traits that mediate the ecological interactions that generate variation in 

fitness, it is impossible to capture all of them in any one study. In a system like ours, where 
682 we employed a novel predator for Drosophila, anti-predator behaviors that were initially 

rare in the progenitor population can rise in frequency, fundamentally changing aspects 
684 of selection on other traits. For instance, if the escape response to direct attacks was the 

primary strategy early in the evolutionary process, but the ability to avoid the predators 
eee developed in later generations, then selection on wing size and shape would potentially 

become far weaker. Study systems like the one used here allow for additional future work 
ess to address these questions in a relatively straightforward manner, which will likely be- 
come increasingly important as we recognize the limitations of measuring relatively small 
ego numbers of traits. 
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Table 1. Survival ability in the base population measured as the number of surviving 
flies in each arena after 24 hours exposure with the predators. Table shows the output 
from the LM funtion in R for a set of models evaluated using Bayesian information 
criteria. Numbers in parenthesis are standard errors of the above estimates. 
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Table 2. Survival ability in the evolved populations measured as the number of 
surviving flies in each arena after 24 hours exposure with the predators. Table shows the 
output from the LM funtion in R for a set of models evaluated using Bayesian 
information criteria. Numbers in parenthesis are standard errors of the above estimates. 
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Table 3. The output from the logistic regression of wing size onto survival in the base 
population for a set of models evaluated using Bayesian information criteria. Numbers 
in parenthesis are standard errors of the above estimates. The logistic regression models 
were used to evaluate statistical significance of the estimated selection differentials. The 
values reported in the manuscript were taken from identical linear regression models 
(output not shown) 
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Table 4. The output from the logistic regression of wing size onto survival in the 
evolved populations for a set of models evaluated using Bayesian information criteria. 
Numbers in parenthesis are standard errors of the above estimates. The logistic 
regression models were used to evaluate statistical significance of the estimated selection 
differentials. The values reported in the manuscript were taken from identical linear 
regression models (output not shown) 
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Table 5. R 2 values from partialleast squares analysis for shape. 
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Figure 1. (A) 1st instar nymph of the Chinese mantid (Tenodera aridifolia sinensis) 
consuming a fruit fly. Note the wing about to drop off. (B) The selective function for 
size estimated by fitting cubic splines (sensu Schluter, 1988) along with estimates for 
linear and quadratic selection. Stars denote significance from logistic regression, but 
estimates are derived from a linear regression of size on relative fitness. Points above the 
function are individuals that survived. Points below the line were captured and eaten. 
Red dots are females and blue dots are males. Error bands are 95% confidence intervals. 
(C) Visualization of the selection differential (S) as measured in the base population. 
Points indicate landmarks and semi-landmarks. The shapes represent the mean shape 
plus lOx S (black line) and minus lOx S (red line). 
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Figure 2. (A) Mean number of survivors in each predation arena after 30 generations 
of experimental evolution. (B) Differences in wing size of the evolved populations after 
30 generations of experimental evolution. Female values shown are 13% larger than 
those for males. Errors are 95% confidence intervals. 
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Figure 3. (A) Shape score by generation for control (B) and predation selection 
regimes. Model adjusted shape score for generation (sensu Drake and Klingenberg, 2008 
see methods) is plotted against generation number, with blue points for males and red 
points for females. Solid regression lines and 95% confidence intervals are for replicate 1, 
and dashed lines and 95% confidence intervals are for replicate 2. (C) Visualization of 
the directions of the evolution of wing shape in the 4 experimental evolution 
populations. The shapes represent the mean plus (black line) and minus (red line) the 
modelled vector of evolutionary change in each case, scaled to 50 generations in 
magnitude. The points represent landmarks and semi-landmarks. Vector correlations 
between these modelled directions of shape evolutions (and their 95% credible intervals) 
are printed between the pairs of populations to which they relate. 
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Figure 4. The selective function for size estimated by fitting cubic splines (sensu 
Schluter, 1988) with replicates pooled for the (A) control (B) and predation populations 
along with estimates for linear and quadratic selection. Points above the function are 
individuals that survived. Points below the line were captured and eaten. Red dots are 
females and blue dots are males. Error bands are 95% confidence intervals. (C) 
Magnitude of the selection differential (S) as measured in the base (B), the control, and 
predation populations. Black points and lines are estimates and bootstrapped 95% 
confidence interval. The grey lines are the 95% confidence intervals from permutation of 
the same data; they represent the null hypothesis that the magnitude of S is random 
relative to survival. 
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Figure 5. Vector correlations between S for wing shape estimated in the base 
population, and the direction of the direction of shape change during experimental 
evolution. The response vector was estimated within each population. Points are vector 
correlation estimates, and lines represent 95% bootstrapped confidence intervals. 
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Figure 6. Estimates of variance for shape calculated as the trace of the covariance 
matrix for female and male flies from the evolved populations. Estimates of total 
variance (grey diamonds) are calculated with dead and surviving flies combined. Error 
bars are 95% bootstrapped confidence intervals. 
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Figure S 1. Proportion of wild-type flies surviving in each arena. Error bars are 95% 
confidence intervals. 
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Figure S 2. A greyscale wing image from a microscope-mounted camera, showing the 
splines fitted by the WingMachine software (left panel). Landmarks and semi-landmarks 
can be extracted from these splines, and superimposed using the CPReader program. 
From this data the wing is modelled by the position of these landmarks and 
semi-landmarks (represented by the open circles in the right panel), and simply 
connecting these coordinates with straight line segments as in the right panel allows 
for an easily interpretable visualisation of the wing veins and margins. 
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Figure S 3. Vector correlations between selection differentials (S) measured in the base 
population (B) and all four experimental evolution populations. Points and lines are 
vector correlations and bootstrapped 95% confidence intervals. 
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Figure S 4. Visualisation of the selection differential (S) as measured at generations 31 
Sz 32 (see methods) in the experimental evolution populations. The shapes represent the 
mean shape plusl5x S (black line) and minus 15x S (red line), and the points represent 
landmarks and semi-landmarks. Vector correlations between S vectors (and their 95% 
credible intervals) are printed between the pairs of populations to which they relate. 



